/*  scftrgbs.cpp  by K.Tsuru  */
// function ID = 3320, 3321, 3322, 3323, 3324 ver. 2.18
/****************************************************************
It evaluates cos(x) and sin(x) using the binary splitting method.

<< prototype in "snfunc.h">>
TrigFuncValues CosSinBS(const SDouble& x, bool getTan = false);
where
    struct TrigFuncValues {
        int quadrant;
        SDouble angle,  cos,sin,tan;
    };
and if set "getTan = true" tangent value is evaluated
 else a request to use/show tangent value yields an error.

CPU time of sin(x) and cos(x) in second
f/4  1,000  5,000* 10,000* 20,000*
SL    0.6    3.1    6.5     13.6
SD    0.6    3.6    6.5     15.6
old   0.2    3.6   12.2     46.3

*SLong version is slightly faster.

Bug fix of CosBS(x), etc. (on June 6, 2015) version 2.21
****************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

static const SLong S_BASE = Lpow10(4*DFIGURES); // consider efficiency of FFT mult.

#define USE_SLC 1   // use SLComplex or not

#if USE_SLC
typedef SLComplex ExpImagBSRType; // SLong version exp(id*pi) 6.30(sec) in 40,000 decimals
#else
typedef SComplex ExpImagBSRType;// SDouble version  7.31 (sec)
#endif
/******************************************************
It returns the value exp(i n/d) as a rational number a/c.
******************************************************/
static SDouble C;
static const char* const funcName = "ExpImagBSR()";
static void ExpImagBSR(const ExpImagBSRType& n, const ExpImagBSRType& d,
        ExpImagBSRType& a, ExpImagBSRType& c) {
    long prec = long(C.EffFig() + C.Hidden()) * DFIGURES;
    double log10x = CeilLog10(n.Imag()) - CeilLog10(d.Real());

    if(log10x > 0) C.SetError(C.SYNTAX_ERR, funcName, 3320);
    long L = upToExpSeries(prec, log10x + 1);
    ExpBSRationalNumber<ExpImagBSRType> exp_ix(L, prec, n, d);  // defined in "sbstempl.h"

    exp_ix.putTogether();
    a = exp_ix.getA();  c = exp_ix.getC();  //  exp_ix.getAC(a, c);
}
/************************************
  cosVal=cos(x), sinVal=sin(x)  3320
  It returns the quadrant of angle x.
*************************************/
int CosSinBS(const SDouble& x, SDouble& cosVal, SDouble& sinVal) {
    SDouble y, r;
    int func = COS_CALC, quadrant;
    int sgn = GetTriCalcMethod(x, y, &func, &quadrant);// |y| < pi/4
    if(sgn == 0){   // y = -1, 0, 1
        x.upToTerm = 0;     cosVal = y;
        func = SIN_CALC;    GetTriCalcMethod(x, sinVal, &func);
        return 0;
    }
// sgn != 0 : cos(x) = sgn*func(y)  (|y|<pi/4 )
    SNBlock <SLFraction> slr;
    int k = SplitSDouble(y, slr, S_BASE);   // very fast

    SComplex rA(1.0, 0.0), rC(1.0, 0.0);
    ExpImagBSRType num, den, a, c;

    for(int i = 1; i < k; i++) {
        num.Set( 0.0, slr(i).num);      den.Set( slr(i).den, 0.0);
        ExpImagBSR(num, den, a, c); // a/c = exp(i*n/d)
#if USE_SLC // SLong Version
        rA *= SLCtoSC(a);   rC *= SLCtoSC(c);
#else       // SDouble Version
        rA *= a;        rC *= c;
#endif
    }
    SComplex result = (rA / rC);

    cosVal = sgn * ( (func == COS_CALC) ? result.Real() : result.Imag() );
    sinVal = (func == COS_CALC) ? result.Imag() : result.Real();
    if (quadrant <= 2 && sinVal.Sign() < 0) sinVal.ChangeSign();
    if (quadrant >= 3 && sinVal.Sign() > 0) sinVal.ChangeSign();
    return quadrant;
}
////// copy constructor and substitution operator //////
TrigFuncValues::TrigFuncValues(const TrigFuncValues& a)
    : quadrant(a.quadrant), angle(a.angle), cos(a.cos), sin(a.sin) {
            if(a.tan.RawSign() != tan.UNDECIDED) tan = a.tan;
}
TrigFuncValues& TrigFuncValues::operator=(const TrigFuncValues& a) {
    if(this != &a) {
        quadrant = a.quadrant;
        angle = a.angle; sin = a.sin; cos = a.cos;
        if(a.tan.RawSign() != tan.UNDECIDED) tan = a.tan;
    }
    return *this;
}
/***************************
main body
x = ip + dp
exp(i*x) = cos(x) + i*sin(x)
****************************/
// initialized by angle = 0, cos=1, sin=0, tan=0
static TrigFuncValues a(0.0, 1.0, 0.0, 0.0);

// cos(x), sin(x) and tan(x) 3321
TrigFuncValues CosSinBS(const SDouble& x, bool getTan) {
    if( a.angle != x ) a.quadrant = CosSinBS(a.angle = x, a.cos, a.sin);
// Here a.angle == x then "a.cos" and "a.sin" ware already calculated.
    if( (a.cos.Sign(3321) == x.ZERO) || !getTan ) a.tan.SizeZero(); // sign = UNDECIDED;
    else a.tan = a.sin * DReciprocal(a.cos);

    return a;
}

SDouble CosBS(const SDouble& x) {// cos(x) 3322
    TrigFuncValues tv = CosSinBS(x, false); // Do not use "a" which is a static member above.
    return tv.cos;
}

SDouble SinBS(const SDouble& x) {// sin(x) 3323
    TrigFuncValues tv = CosSinBS(x, false);
    return tv.sin;
}

SDouble TanBS(const SDouble& x) {// tan(x) 3324
    TrigFuncValues tv = CosSinBS(x, true);
    return tv.tan;
}